The L-shaped selection algorithm for multitrait genomic selection

Abstract Selecting for multiple traits as opposed to a single trait has become increasingly important in genomic selection. As one of the most popular approaches to multitrait genomic selection, index selection uses a weighted average of all traits as a single breeding objective. Although intuitive and effective, index selection is not only numerically sensitive but also structurally incapable of finding certain optimal breeding parents. This paper proposes a new selection method for multitrait genomic selection, the L-shaped selection, which addresses the limitations of index selection by normalizing the trait values and using an L-shaped objective function to find optimal breeding parents. This algorithm has been proven to be able to find any Pareto optimal solution with appropriate weights. Two performance metrics have also been defined to quantify multitrait genomic selection algorithms with respect to their ability to accelerate genetic gain and preserve genetic diversity. Computational experiments were conducted to demonstrate the improved performance of L-shaped selection over-index selection.


Introduction
The effectiveness of genomic selection (GS) in accelerating genetic gain in plant and animal breeding programs (Jannink et al. 2010;Meuwissen and Goddard 2010;Desta and Ortiz 2014;Rutkoski et al. 2016) has motivated breeders to apply the technique for multiple traits, including yield, quality, and tolerance to biotic and abiotic stresses (Jia and Jannink 2012). Breeders are more likely to invest in GS programs that are capable of improving multiple traits of a crop rather than a single trait throughout the generations (Lynch and Walsh 1998;Bernardo 2002).
Existing approaches for multitrait genomic selection (MTGS) include tandem selection, independent culling selection, and index selection. Tandem selection treats MTGS as the aggregation of multiple single-trait GS programs and selects for the traits sequentially (Burgess and West 1993). Independent culling selection sets a minimum threshold (i.e. culling levels) for each trait and only selects individuals that exceed the culling levels for all traits (Lorenzana and Bernardo 2009). Index selection converts MTGS to a single-trait GS by using a linear combination of individual traits weighted by their importance as the breeding objective (Hazel and Lush 1942;Hazel 1943;Williams 1962). Recently, Moeinizade et al. (2020) proposed a new algorithm for MTGS that maximizes one trait subject to the constraints that another trait falls within a desirable range.
Index selection has heretofore been a commonly used approach to MTGS due to its capability to select for multiple traits simultaneously and its flexibility to assign different weights according to the relative importance of the traits. In contrast, tandem selection can only select for one trait at a time, and independent culling selection may eliminate an otherwise highperforming individual due to its minor shortcoming in one trait (Lorenzana and Bernardo 2009). Index selection overcomes such limitations by taking an importance-weighted linear combination of all traits, giving breeders a wide range of trade-off options among the traits to choose from.
However, index selection also suffers from its own limitations, i.e. its numerical sensitivity and inability to find certain optimal selections. In this study, we propose a new approach to MTGS, which uses an L-shaped objective function (as opposed to the linear objective function used in index selection) to select optimal breeding parents that strike a balance among multiple traits with respect to their relative importance. This algorithm not only overcomes the 2 limitations of index selection but also demonstrates superior performance with respect to both accelerating genetic gain and preserving genetic diversity.
The rest of the paper is organized as follows. In Materials and methods, we formally formulate MTGS as a multiobjective optimization problem, introduce the L-shaped selection algorithm, and present the mathematical properties that allow it to overcome the limitations of index selection. Moreover, we define 2 metrics for assessing the performances of MTGS algorithms in terms of accelerating genetic gain and preserving genetic diversity.
In Computational experiments, we describe the computational experiments that we conducted to compare index selection and

Materials and methods
Consider a breeding project that starts with an initial population of plant or animal individuals. A number of crosses are made in each generation to produce a new population of progeny in the next generation until a predefined deadline for the project. Suppose there are multiple traits that the breeders aim to improve through the breeding process. Under the following 3 simplifying assumptions, the focus of our study is to select the right individuals to make the right crosses in order to optimize all traits at the end of the breeding project.
Assumption 1 There are adequate and reliable genetic data, including genotype and recombination frequencies.
Assumption 2 All traits are largely determined by additive effects with negligible dominance or epistatic effects.
Assumption 3 Allele effects for all traits have been estimated sufficiently accurately and constant over the breeding process.

Problem definition
The objective of MTGS is to select a subset of breeding parents from a group of candidate individuals based on their genotype and estimated allele effects in order to maximize genetic gains with respect to multiple traits over a number of breeding generations. The following nomenclature will be used in this paper.
I set of candidate individuals of plants or animals for selection J set of loci K set of traits G i;j genotype of individual i 2 I at locus j 2 J b j;k effect of allele j 2 J on trait k 2 K v i;k genetic value of individual i 2 I on trait k 2 K: v i;k ¼ P j G i;j b j;k w k weight parameter that indicates the relative importance of trait k 2 K x i binary variable indicating whether individual i 2 I is se- number of breeding parents to be selected Without loss of generality, we assume that maximization (rather than minimization) is the direction of improvements for all traits. For a trait k that needs to be minimized, we can replace v i;k with Àv i;k for all i 2 I in a maximization model, which is equivalent to minimizing trait k. For traits whose values need to be contained within a desirable range, we can maximize the percentage of individuals in a population whose trait values fall within such a range.
With the above definitions and assumptions, the MTGS can be formulated as the following multiobjective optimization model: x i 2 f0; 1g 8i: Here, the objective function (1) is the maximization of all traits of selected individuals. Constraint (2) requires that exactly S individuals be selected. Constraint (3) defines x i as a binary variable for all i 2 I.
In multiobjective optimization, a feasible solution is called Pareto optimal if it is not dominated by any other feasible solution. Solutionx dominatesx ifx is no worse thanx in any trait and better in at least one: Moreover, the ultimate goal of solving a multiobjective optimization problem is to find not a single Pareto optimal solution but all Pareto optimal solutions that represent the range of possible trade-offs among different traits, referred to as the Pareto frontier. Solutions on the Pareto frontier are not equivalent; rather, they offer nondominated alternatives for decision-makers to choose from, with each alternative being optimal under certain criteria of trade-off among the multiple objectives.

Index selection
As a widely used selection approach for MTGS, index selection solves a single objective optimization model that maximizes the weighted average of all traits: Here, the objective function (4) is the weighted average genetic value of all traits. When different weight parameters are used, index selection essentially searches for the convex efficient frontier, which is the subset of Pareto optimal solutions that are on the convex hull of the feasible region.
The following proposition shows that, when strictly positive weights are used for all traits, the optimal solution to index selection must be Pareto optimal to MTGS.

Two limitations of index selection
Despite the effectiveness of index selection in finding Pareto optimal solutions to MTGS, it suffers from 2 major limitations. First, numerical solutions to (4)-(6) may be sensitive to the units being used for the genetic effects of different traits. For example, when we try to maximize both plant height and grain yield with their respective weights, different solutions may result from simply changing the units used to measure the 2 traits from inches and bushels per acre to meters and tonnes per hectare. As a result, it is challenging to determine or interpret the numerical values of the weights for different traits. The second limitation is the inability to discover all Pareto optimal solutions by using different values of weight parameters w k due to the convexity of the objective function (4) and the potential nonconvexity of the set of Pareto optimal solutions. To illustrate this point, consider the following example.
Example 1. Suppose 2 individuals are to be selected from 4 to make a cross, then there are 6 candidate crosses. The genetic values of the 4 individuals and 6 crosses for 2 important traits are summarized in Tables 1 and 2. All 6 crosses are Pareto optimal since none of them is dominated by another. However, as shown in Fig. 1, the index selection model (4)-(6) can only find 3 crosses (c 1 , c 5 , and c 6 ) that are on the efficient frontier of the 6 crosses in the v 1 -v 2 space. The other 3 will never be selected no matter what non-negative weights w 1 and w 2 are used because they are dominated by the line segment between c 1 and c 5 .

L-shaped selection
The proposed L-shaped selection for MTGS can be formulated as the following optimization model: is the normalized genetic value that falls within the range of (0, 1), where v k and v k are the lower and upper bounds of trait k, which can be obtained, respectively, as v k ¼ min i ðv i;k À Þ; 8i 2Ĩ and v k ¼ max i ðv i;k þ Þ; 8i 2Ĩ for a large set of individuals I ; here is a small positive value to ensure that the normalized genetic value falls within (0, 1) and not on the boundaries.
This new formulation was designed to address the 2 limitations discussed in Two limitations of index selection. The normalized genetic effectṽ i;k is independent of the measurement units being used. The application of this well-known normalization technique in MTGS allows the model to strike a balance among multiple traits in more meaningful terms. After normalization, the numerical values of the weights can be interpreted as the relative importance for each percentage point of improvement in the respective traits. For example, if w 1 ¼ 0:8 and w 2 ¼ 0:2, then a 1% improvement of trait 1 is 4 times as important as a 1% improvement of trait 2, regardless of the measurement units used for the 2 traits. This interpretation of weights allows breeders to make the most appropriate trade-off among different traits. In certain cases, relative monetary values of the traits can be used as the weights; in other cases, such as when the prices of certain traits are volatile or otherwise hard to determine, the weights can also reflect subject emphases that breeders place on different traits.
The following proposition demonstrates how model (7)-(9) addresses the second limitation.
Proof. We claim that w k :¼ P ixiṽi;k > 0; 8k 2 K are such thatx is optimal to (7)-(9). Suppose not, then there exists a solutionx such that min k P ixiṽi;k w k > min k P ixiṽi;k w k ¼ 1: Fig. 1. The left subfigure shows the 6 possible crosses in the v 1 -v 2 space and the efficient frontier that index selection is able to find. The right subfigure shows that only c 1 , c 5 , and c 6 could be found to be Pareto optimal using index selection with different weights w 1 and w 2 .
ixi v i;k ; 8k 2 K: This contradicts the assumption thatx is Pareto optimal to (1)-(3). Therefore,x must be optimal to (7)-(9) with w k :¼ P ixiṽi;k ! 0; 8k 2 K. h Proposition 2 is a guarantee that L-shaped selection is able to find any Pareto optimal solution to MTGS with appropriate weight parameters. To illustrate this desirable property, which index selection does not have, we solved Example 1 using Lshaped selection, and results are shown in Fig. 2. The left subfigure illustrates that, in the case of k ¼ 2, model (7)-(9) is trying to slide an L-shaped objective function (hence the name) along the direction from the origin toward (w 1 , w 2 ) to the maximal extent while touching at least one solution with the Lshaped curve. The right subfigure shows the different Pareto optimal solutions that can be found using different combinations of weight parameters w 1 and w 2 .
The significance of Proposition 2 is the ability of L-shaped selection to offer a diverse set of alternative breeding solutions on the Pareto frontier that empowers breeders in 3 important ways: (1) understand the set of available trade-offs among multiple traits, (2) realize the outcomes as a result of a wide range of weights being used, and (3) select the most appropriate weights and breeding parents according to the specific breeding objectives and constraints.

Performance measures of MTGS algorithms
In this section, we define 2 measures, namely Pareto optimality gap and diversity, to evaluate the performance of MTGS algorithms. The motivation is to assess the capability of an algorithm to produce progeny through the breeding process that is not only Pareto optimal but also representative of diverse trade-offs among different traits.
Suppose an MTGS algorithm was used in a number of breeding projects, each with different sets of weight parameters. Let I 0 denotes the set of individuals produced in the final generation of all breeding projects combined, and let I 1 denotes a superset of I 0 , possibly also including individuals produced from all other competing algorithms. For any set of individualsÎ , we define PðÎ Þ as the Pareto optimal subset ofÎ : The 2 performance measures are defined as follows: • Pareto optimality gap of PðI 0 Þ against PðI 1 Þ is defined as which measures the extent to which set PðI 0 Þ is dominated by PðI 1 Þ. The term ðv i 0 ;k À v i;k Þ þ :¼ maxfv i 0 ;k À v i;k ; 0g detects any positive gap between individuals i 0 2 PðI 1 Þ and i 2 PðI 0 Þ on trait k. Then, the term min k2K w i;k ðv i 0 ;k À v i;k Þ þ identifies the smallest weighted gap across all traits in order to determine the extent to which individual i is dominated by i 0 . Next, the term min i 0 2PðI 1 Þ min k2K w i;k ðv i 0 ;k À v i;k Þ þ identifies the individual i 0 that dominates i by the least amount. This term is the Pareto optimality gap between individual i and PðI 1 Þ, and the summation of which over all individuals in PðI 0 Þ gives the Pareto optimality gap of the set PðI 0 Þ against PðI 1 Þ. The Pareto optimality gap can be considered as the minimal amount of trait improvements in I 0 necessary to break the dominance of PðI 1 Þ over any individual in PðI 0 Þ. A noteworthy observation here is that the Pareto optimality gap between i 2 PðI 0 Þ and i 0 2 PðI 1 Þ is zero if v i 0 ;k 0 ¼ v i;k 0 for one trait k 0 and v i 0 ;k > v i;k for all others k 2 P k 0 . This may be counter-intuitive but is also defensible because, although i is dominated by i 0 , it takes an arbitrarily small positive improvement in individual i on trait k 0 to break the dominance. Fig. 2. The left subfigure shows the 6 candidate solutions in the v 1 -v 2 space and how L-shaped selection uses an L-shaped objective function to search for Pareto optimal solutions. As an example, c 3 is found to be optimal with equal weights on the 2 traits, since it allows the L-shaped objective function to slide the furthest away from the origin toward the direction ðv 1 ¼ w 1 ¼ 0:5; v 2 ¼ w 2 ¼ 0:5Þ. Moreover, c 3 is optimal for all weights inside the shaded (unbounded) triangle, the 2 edges of which cross the vertices of 2 L-shaped objective functions with one crossing c 2 and c 3 and the other crossing c 3 and c 4 . The right subfigure shows that all 6 candidate solutions can be found to be Pareto optimal using L-shaped selection with different weights w 1 and w 2 . The 6 candidate solutions are also plotted in the right subfigure in the w 1 -w 2 space to illustrate how different regions of weights are determined.
• Diversity of PðI 0 Þ is defined as which measures the weighted average Euclidean distance of all Pareto optimal solutions within set PðI 0 Þ. An ideal MTGS algorithm should be able to produce progeny with a small Pareto optimality gap (not being dominated by progeny produced from competing algorithms) and a large diversity (offering different trade-off options in traits).

Computational experiments
We compared the performances of index selection and L-shaped selection with computational experiments using a maize data set considering 2 traits: plant height and ear diameter.

Data sets
We used 200 maize inbred lines of 369 shoot apical meristem population distributed across the 10 chromosomes (ISU 2020). We extracted 1,000 single nucleotide polymorphisms (SNPs) out of the total of 1.4 million SNPs that were collected using genomewide association study used by Leiboff et al. (2015), merged with additional SNPs genotyped using tGBS (Schnable et al. 2013), and those which were phased using Beagle (Browning and Browning 2009). Recombination rates in this population were estimated using the genetic map developed from the maize nested association mapping population. Genetic effects of the 2 traits, plant height and ear diameter, were extracted from Bernardo and Yu (2007). Both traits are affected by 1,000 loci with pleiotropic effects of varying magnitudes. The genetic correlation of the 2 traits is 0.375. Like index selection, L-shaped selection also applies to both correlated and uncorrelated traits. A similar data set was used in Amini et al. (2021), which also showed figures of summary statistics of the data.

Breeding process
In each simulation of the breeding process, 200 individuals were randomly selected from the 369 inbred lines to form an initial population. In each of the subsequent generations, 2 individuals were selected using either index selection or L-shaped selection to produce 200 progeny in the next generation. The genetic values v i;k of all 200 individuals i for the 2 traits k in the fifth generation were used for performance analyses. Nine groups of this breeding process were simulated, each for a different set of weight parameters w 1 2 f0:1; 0:2; . . . ; 0:9g and w 2 ¼ 1 À w 1 with 10 independent repetitions.

Results and discussions
We compared the performances of index selection and L-shaped selection with respect to normalized plant heights and ear diameters of progeny in the final generation (T ¼ 5), with both traits to be maximized. Figure 3 shows the aggregated results of 90 experiments (9 weights by 10 repetitions) using 2 algorithms; only those progeny that were Pareto optimal within each experiment were plotted. The x-axis and y-axis represent the normalized plant height and normalized ear diameter, respectively. It can be seen that L-shaped selection resulted in better-performing progeny in terms of both genetic gain and diversity. In fact, the same number of progeny was produced using the 2 selection methods. The reason that the left subfigure had much fewer points than the right one was because, compared with L-shaped selection, much more progeny from index selection were either dominated by or overlapping with other progeny in the same generation.
The optimal Pareto frontier of Fig. 3 is shown in Fig. 4. It can be seen that, for all different sets of weight parameters w 1 2 f0:1; 0:2; . . . ; 0:9g and w 2 ¼ 1 À w 1 , L-shaped selection outperforms index selection, since the optimal Pareto frontier for Lshaped selection not only dominates that for index selection but also is more diverse. Table 3 shows the Pareto optimality gap and diversity of all progeny in the final generation. The Pareto optimality gap assesses the capability of a selection method in resulting in Fig. 3. Performance of progeny in the final generation for index selection (left) and L-shaped selection (right).
Pareto optimal progeny; the lower the gap, the more Pareto optimal the progeny. We assessed the Pareto optimality gap of index selection and L-shaped selection against the combined sets of progeny from the 2 approaches. As can be seen in Table 3, the Pareto optimality gap for L-shaped selection is zero, meaning that none of its Pareto optimal progeny was dominated by those from index selection. In contrast, index selection had a Pareto optimality gap of 0.2916, meaning that the Pareto optimal progeny from index selection was on average dominated by 0.2916 by those from L-shaped selection. The table also shows that progeny in the final generation produced using L-shaped selection were more diverse than those using index selection, in terms of representing different trade-offs between the 2 traits. These observations were consistent with results from Fig. 3.

Conclusion
We presented the L-shaped selection for MTGS and demonstrated its improvements over the commonly used index selection in the capability to produce elite progeny with better genetic traits and higher diversity. Motivated by 2 limitations of index selection, the new approach made 2 major contributions to MTGS. First, L-shaped selection is robust against different measurement units used for multiple traits. Second, L-shaped selection guarantees to find all Pareto optimal solutions with appropriate weight parameters. Moreover, we introduced 2 metrics to quantify the capability of MTGS algorithms in accelerating genetic gain and preserving genetic diversity.
Besides theoretical contributions, L-shaped selection also outperformed index selection in computational experiments using a maize data set. The results demonstrate that L-shaped selection outperforms index selection in producing better-performing and more diverse population, considering plant height and ear diameter as 2 traits. Supported by the theoretical results, L-shaped selection is also applicable to other breeding populations with available genotypic and phenotypic data sets.
This study is not without limitations that could be addressed in follow-up research studies. Firstly, one can integrate the objective function in the L-shaped selection formulation in more sophisticated algorithms for GS (De Beukelaer et al. 2017;Goiffon et al. 2017;Gorjanc et al. 2018;Moeinizade et al. 2020) rather than a straightforward truncation selection as used in this study. Secondly, although we have considered only 2 traits in our study in order to better visualize the results in 2D figures, the method applies to an arbitrarily large number of traits. A comprehensive case study with appropriate data from a large number of traits would be of interest for a future research project. Thirdly, Lshaped selection was designed to improve the simple nonpenalized index selection. In contrast, the multitrait look-ahead selection (MT-LAS) algorithm (Moeinizade et al. 2020) was designed to specifically address the case where one trait is to be maximized and the other restricted between a lower bound and an upper bound, which was comparable with the penalized index selection (Lopez-Cruz et al. 2020). However, MT-LAS is much more computationally intensive and not compatible with other MTGS methods, whereas L-shaped selection can be combined with other MTGS methods by replacing their objective functions with the Lshaped one. Finally, it would be more realistic to consider dominance, epistatic, and environmental effects in calculating genetic trait values of individuals (Amini et al. 2021), however, only additive effects have been considered in this study.

Data availability
The datasets used in the computational experiments were derived from sources in the public domain as described in Data sets. Computer codes used for the computational experiment were uploaded to GitHub (2022).